Soumya D. Sanyal

In this post, I'm going to try and see if I can elicit any relationships between the explanatory features and the explained features (expenditure).

This part of the project serves to aid in getting to actually know the data, and in hypothesis making, hypothesis refinement, and model selection; it essentially lays the foundation for the "story" that the project will tell.

I find this stage of the modeling process to be one of the most interesting. It's practically impossible to know if your priors about the problem you're looking at will influence the modeling process - I think they almost certainly will. At the same time, there's dramatic tension - will my hunch about the problem be supported by the data or will the data show me that it's completely false? Will another insight rise to take it's place? Will there be any story at all? All of these outcomes are entirely possible at this stage, and it's quite nerve wracking.

In any case, I can do my best by examining the data for as much supporting (or refuting) evidence in as objective a manner as possible.

Here are some of the questions I'm hoping to answer at the end of this post:

1. What explanatory features are predictive (if any)?
2. What explains the statistical variation in the dataset? What does the randomness look like? Is there a pattern? Surely it's not all just white noise!
3. Which features are mutually correlated, and which are mutually uncorrelated? Do any features combine to create stronger relationships?

Next, we've accumulated a lot of boilerplate code that I need to load to start up my modeling session at this point. Rather than subject you to a lot of scrolling, I'm simply going to save it to a python module called blogging.py, and import it in one step. If you'd like to take a look at the environment I'm importing to build the analysis in this post, it's here.

The next step in my analysis is to explore the data a bit more. I want to see if I can get any more insight into what the underlying data generating process is; the more I can rely on a clean look at the data, the better I'll be able to dispel my own incorrect assumptions about what is or should be happening.

One problem I do have to contend with in working with this dataset is that it comes as weighted survey data. Working with weighted data can be quite tricky. In the very first instance, since the unweighted data does not (by definition) accurately represent the population of interest, it's impossible to draw valid inferences from the standard exploratory (or even regression) procedure without accounting for the weights.

When constructing summary data, as in a histogram, or constructing a linear function of the data, this is not much of a problem. Any linear function $f(x)$ of the data is just as easy to compute as $f(Wx)$, where $w$ is a weight matrix, and variances will just be quadratic in $W$ as expected. However, there are myriad other challenges in working with weighted data; there's a relatively recent discussion of some of these challenges in this paper by Andrew Gelman.

For my immediate purposes - my Insight project - the challenge was to come up with some rough and ready ways to accommodate these weights. While the histogram methods in the Python ecosystem have options for specifying item level weights, not all EDA function calls do. In particular, scatterplots don't (as far as I know). Most model estimation calls that I looked into before starting this project did have some accommodations for item weights, so I figured that if I could get past the EDA stage, I might be fine.

I needed to find a way of getting a sense of the pairwise correlations in the dataset while accounting for the weights. In the end, I found that matplotlib hexbin plots had some support for item weights, and I figured that adding the weighted sum of points to hexbin cells would do well enough without my having to write my own visualization functions. In the end, this is what I constructed:

Let's try this out on a an example pair of features.

In this graph, greater (weighted) mass in any one hexbin results in a higher color count (according to the colorbar scale on the right). More tightly packed data results in a greater plot density.

There looks to be a general tendency for older people to spend more on office based services in this picture, but my sense of the trend is obscured by the outliers above 100k. Of these outliers, most are in the 70-80 year range, with one below 20 years. There's not enough data for me to make a robust case about whether age correlates well with the incidence of really high expenditure, so I'm going to try looking at this again after cutting off the outliers.

This is helpful. It looks like a lot of weighted mass (the aqua region at the bottom) is concentrated near zero expenditure. However, it highlights the fact that expenditure trends higher with age (as does variance in expenditure).

It looks like this is a fairly useful visualization technique as far as quick and dirty goes.

Here are the other groups of correlations I want to inspect:

1. Income
2. Utilization
3. Health conditions
4. Duration without insurance
5. Age of diagnosis
6. BMI
7. Demographic information
8. Expenditures.

I also want to look at the four expenditure categories at the same time. Here's a quick method to help me do that:

It looks like the same pattern holds for the other expenditure categories as for office expenditure when it comes to age.

Looking at income:

Here you get the result that while the level of spending moves away from zero as income rises, the variance of expenditure in all four categories declines with income. An important question is whether the reduced variation in spending at higher incomes is simply a result of having fewer people at those income levels, or if there is a systematic effect (more effective health management, better bargaining power).

Looking at how utilization correlates with expenditure:

I get the following out of these plots:

1. Office based utilization and expenditure are positively correlated. Office based utilization is negatively correlated with all other expenditure.
2. Outpatient utilization doesn't strongly correlate with outpatient expenditure. Outpatient utilization correlates somewhat negatively with other categories of expenditure.
3. ER visits correlate negatively with other expenditures, and not particularly at all with ER expenditure.
4. Inpatient visits correspond positively with inpatient utilization, and negatively with all other expenditure.

Next, I'll look at chronic health diagnoses. There are over 40 in my list, so it'll be useful to break them into groups.

Before we take a look at these plots, I want to quickly mention a very counterintuitive finding that I came across when I looked at these. For many of these chronic conditions, the subpopulation that does not have a positive diagnosis has a much heavier right tail than the subpopulation with the chronic condition.

On the face of it, it looked like people overall were at greater risk of spending more if they did not have a chronic illness than if they did. The paradox disappears when you zoom in on these plots around the first $\$10,000$ spent; it shows that many more people without a chronic condition are clustered at $\$0$ spent, relatively speaking. It's interesting to conjecture what may cause the thinner tails for chronic sufferers - one explanation may just be that the larger subpopulation is naturally accompanied by a larger variance.

The plots below shows the scatterplots to this limited scale for emphasis.

Conclusions:

1. A cancer diagnosis overall drives expenditure for office based services and inpatient services, and has less of a clear effect on outpatient spending and ER spending. For these latter two categories, there is greater variance in spending for the subpopulation without a positive cancer diagnosis.
2. There is very little mass in the graphs measuring spending against the various subdiagnoses of cancer. It seems likely that there will be a corresponding lack of statistical significance of these features in any analysis I do, and it may be well to leave those out and just use the overall cancer diagnosis feature.

Conclusions:

1. Diagnoses of emphysema, bronchitis, and asthma drive spending on office services and inpatient services. They have less of an effect on outpatient and ER services. The right tails in these latter categories is heavier in the negative group.
2. The other features do not seem to be as predictive of increased spending.

Conclusions:

1. Diagnoses of high blood pressure, high cholesterol, coronary heart disease, angina, heart attack, other heart disease, stroke, and diabetes drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.
2. A multiple diagnosis of high blood pressure is probably not very predictive.

Conclusions:

1. Diagnoses of joint pain and arthritis drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.

Conclusions:

1. Diagnoses of ADD/ADHD, limitations in physical functioning, and incidence of pregnancy drive spending on office services and inpatient services, while having a limited effect on the other two categories of spending, in which they are associated with thinner right tails.

Conclusion:

1. It doesn't look like spending is correlated with time elapsed without health insurance.

Conclusion:

1. There doesn't seem to be a strong correlation between the age of diagnosis for any of these chronic conditions and annual spending on any category of healthcare.

Conclusions:

1. No features except insurance status generate significant differences in outpatient expenditure.
2. ER spending is affected by race, marital status, employment status, self-employment status, holding multiple jobs, union status, whether you are offered health insurance by your employer, your overall insurance status (notably whether you have private insurance or a public option), whether your household uses an FSA, where you live (notably, there is more zero clustering for ER spending in the South, and less in the Northeast), and educational attainment.
3. Both office and inpatient spending are affected by race, marital status, employment status, having more than one job, being self employed, union status, whether health insurance is offered through employment, whether public or private options for health insurance are used, whether FSA's are used, insurance status, census region, and educational attainment.

Finally, a look at whether body mass index affects spending:

From looking at this, there doesn't seem to be a strong correlation between spending and BMI in any category. I find this somewhat counterintuitive, since obesity is supposed to be correlated with poor cardiovascular health. I'm inclined to leave this feature in to see if there are any interesting interaction effects that may turn up in the modeling process.

At this point we have a better view of the data, and a sense of what features we expect to be predictive. This looks like a good foundation for a pretty rich modeling exercise.

In the next post, I'll write about the model selection process.

In this post, I'm going to look at the explanatory features in the dataset.

First, boilerplate:

Let's get our dataset, variables and functions from before.

Next I'm going to prototype a function that will spit out the empirical distributions of each of the explanatory variables.

I'm going to have to treat categorical variables differently to continuous variables.

One rough and ready way to break out the categorical variables is to check how many levels that feature has in the dataset - if there are relatively few (say below a threshold of 10 or 20), I'll assume the feature is categorical, and otherwise I'll assume it's continuous. To deal with some of the edge cases, it's easy to note that anything involving an age, or a # of visits should be continuous.

Let's see if this works.

Let's see if this gave us what we wanted.

Pretty good!

Next I'll prototype a function to produce bar plots for the categorical variables. Since the data points are weighted, I need to make sure I scale each point out correctly before I aggregate under each label as well.

Here's the function:

The categorical features above split up nicely into two groups that are interesting. The first consists of health status features, and the second consists of demographic and socioeconomic features. I'll group them accordingly:

Here are a few observations about these health status features:

1. High blood pressure, diabetes, high cholesterol, joint pain and arthritis are relatively common in the population.
2. Heart disease (angina, heart attack, etc) and stroke are less prevalent, but still substantial.
3. Ditto for lung disease (emphysema, bronchitis, asthma).
4. While a cancer diagnosis is relatively common, the individual cancer diagnoses have very low mass in the "yes" bucket. This will likely pose a problem for estimating models, unless we either upsample, or otherwise penalize the other responses. Including these features in our model may be informative if different cancers affect expenditure in different ways, but we'll have to balance this against the lack of reliability of parameter estimates that we'll get from these features. Unless we have a good reason to believe that different cancers cost very different amounts to treat, we may be best off using the aggregate cancer diagnosis feature, and dropping the individual diagnosis features.
5. Practically no-one uses a peak flow meter.
6. There's relatively incidence of ADD/ADHD in the population.
7. Relatively few people are pregnant at the time of the survey.
8. General limitation in physical functioning is fairly prevalent.

Next let's look at the demographic and socio-economic features.

About the demographic/socio-economic features:

1. The male/female subpopulations are evenly balanced in the population, with slightly more women than men.
2. The race/ethnicity distribution is consistent with other surveys.
3. The distribution estimates that about 70m poeple are under 16, a similar number has never been married, about 130m people are married.
4. About 150m people were employed when the survey was conducted. About 100m people are not employed (not in the labor force, or unemployed). This is consistent with reported data elsewhere. About 20m people are self employed, and about 10m people hold more than one job.
5. Remarkably few people - less than 50m - have a choice of health plans at work. I'm interested in whether people who have a choice of health plans tend to pay less for healthcare on average, assuming there is any effect coming from encouraging people to think about their purchasing choices on subsequent utilization or on prices paid/negotiated.
6. The modal entry for industry group and occupational group (by far) is inapplicable. This is a constructed feature in the dataset, so I assume this wasn't directly asked of respondents. It's hard to say what this might mean, or if it'll be useful to keep in our model.
7. About 10m people are in a union. It would also be interesting to know whether being part of a union affects utilization or expenditure in any way.
8. About 60m people did not hold health insurance at the time. These numbers are consistent with pre-ACA uninsured levels.
9. About 80m people were offered health insurance through their employers, and 40m were not.
10. About 180m people under 65 had private insurance coverage during the survey, about 50m people under 65 had only public insurance options, and about 50m people under 65 were uninsured. There were about 20m people over 65 holding a mix of private insurance and medicare, and a similar number with only medicare coverage.
11. About 30m people lived in a household where someone had a flexible spending account. This is about $10\%$ of the population, but FSAs (especially coupled with high-deductible plans) are an interesting demand-side innovation and it would be interesting to understand if this plays a role in determining expenditure.
12. A negligible number of people (relatively speaking) had coverage under a state insurance plan.
13. The distribution of people over census regions is approximately as follows: about 120m in the south, about 75m people in the west, about 65m people in the midwest, and about 60m people in the northeast.
14. The distibution over educational attainment is as follows: about 70m people have an associate's degree past high school, about 65m have a GED or high school diploma, a little under 60m have an 8th grade qualification or less, about 30m people stopped high school between 9th and 12th grade and do not have a high school diploma, a little over 40m people have a 4-year degree, about 25m people have a masters, doctoral or professional degree, and about 25m people fall into the inapplicable bucket.

Now that we've had a look at the categorical features, let's take a look at the continuous features.

There are broadly speaking, continuous variables in the following categories: ages, income levels, time since some event, number of occurences of an event, and total expenditures.

Let's group these variables:

Let's look at the ages first.

The documentation only mentions that the value $-1$ is recorded for age when there is no valid data collected and no imputation can reasonably be made. In the following graph I'll censor this level, since there's not really a good way to interpret this that I can think of. I will probably also leave it out when building the model.

The first thing to notice is the huge mass at 85. There aren't actually over 6m people in the country of age exactly 85; the data is top-coded at this value to preserve respondent privacy. This means I won't be able to build a model that discriminates on age past 85, and I'll just have to lump everyone older than that into the same bucket. I'll hopefully be able to get a reasonable cross section of that subpopulation using other features.

Next let's look at the age of diagnosis features. Since most people don't have chronic illnesses of these types, the mass at -1 dominates the distribution. Accordingly, I'm going to censor those points in these plots, so each of these distributions are conditional on having that chronic condition.

Here's a modified function from the last post for producing histograms with weighted data.

1. The "lifestyle" diseases: heart attack, angina, high blood pressure, coronary heart disease, stroke, emphysema, high cholesterol, and diabetes all follow very similar distributions, which are all not-quite normal, with modes in the 50-60 range, and with some right skew.
2. It will be interesting to look at correlations between these lifestyle disease features. I don't know a priori if this should be high or low, but I would imagine that having multiple comorbidities would be predictive of higher expenditure on healthcare.
3. Arthritis follows a similar distribution to the lifestyle diseases.
4. Asthma is very heavily skewed left, and is most diagnosed at very young ages, with the diagnosis rate trailing off as the population gets older.
5. The age of diagnosis for ADD/ADHD is very tightly clustered in the 0-20 range, with a bell shaped distribution, non-existent tails, and some skew to the left.
6. A generally interesting hypothesis about the relationship between expenditure and the age of diagnosis of any of these diseases may be that the earlier the diagnosis, the better the patient is at making lifestyle changes that can mitigate expenditure levels. It might be interesting to see if this holds up in the modeling.

Let's look at BMI.

The BMI (body mass index) is computed by dividing a person's weight in pounds by their squared height in square inches and multiplying by 703.

Here's how to interpret the adult body mass index according to the CDC guidelines for interpretation.

The distribution of BMI in the population is estimated by this histogram to be bell-shaped, with a heavier right tail, and centered around 25 or so.

Next, on to the income variable:

What's really weird about this is that I see responses below 0; far below 0!

Let's take a closer look.

Thankfully there are not many data points here. An important question is how to interpret an income of $-258220$. The smaller values for negative incomes are also interesting - are these reported after something like adjusting for public assistance by the respondents?

Let's take a look at the survey documentation.

FAMINC13 contains total family income for each person’s CPS family. Family income was derived by constructing person-level total income comprising annual earnings from wages, salaries, bonuses, tips, commissions; business and farm gains and losses; unemployment and workers’ compensation; interest and dividends; alimony, child support, and other private cash transfers; private pensions, IRA withdrawals, social security, and veterans payments; supplemental security income and cash welfare payments from public assistance, Temporary Assistance for Needy Families, and related programs; gains or losses from estates, trusts, partnerships, S corporations, rent, and royalties; and a small amount of “other” income. Person-level income excluded tax refunds and capital gains. Person-level income totals were then summed over family members, as defined by CPSFAMID, to yield CPS family-level total income (FAMINC13). 

Looking at the expenditure features, there's nothing particularly outlandish about these records apart from the incomes. Given this, I think I'm inclined to go by: "when in doubt, throw it out" when it comes to these data points.

Let's truncate at 0 and take a look at the shape of the distribution.

The distribution is unimodal (mode around $\$50,000$), and with most of the mass concentrated between $\$0$ and $\$200,000$, and a thin (but very well off!) tail past that point.

Let's move on to the durations variables.

Since the majority of respondents did carry insurance in the previous year, I'll drop the points at -1 to get a sense of the distribution conditional on not being insured.

This is quite interesting. Here's what the documentation has to say about this feature:

For persons who were covered by health insurance on January 1st, it was ascertained if they were ever without health insurance in the previous year (NOINSBEF). The number of weeks/months without health insurance was also ascertained (NOINSTM, NOINUNIT). 

I'm interested to know if this feature might be predictive of expenditure. It's hard to say - going a long time without insurance might be reflective of being unemployed, or if the respondent did have the means to buy insurance, may be a proxy for their (lack of) demand for risk mitigation. If either of these correlations hold, this might say something about the effect of either bargaining position, or risk tolerance, on propensity to spend.

Let's take a look at this at only positive values.

Next, the utilization features:

It seems clear that most people use relatively little healthcare. Let's confirm if a lot of these utilizations are at 0.

That's definitely the case. Let's see what the distributions look like conditional on positive utilization.

These distributions all look very similar; the modal positive utilization level is one visit, and utilization rapidly decays, although there are very heavy tails for office based services and inpatient services (we've noticed this already in the expenditure features for these categories).

The main difference between these categories is the degree of use. For example, the number of people in the modal office visit level is over 50m, for outpatient and ER is about 30m, and for inpatient nights is about 5m.

Similar pictures in this case.

We've already looked at the patterns in expenditure, so this takes us to the end of the continuous variables we want to look at.

The next thing I'd like to do is to take a closer look at the features to check for bivariate relationships. Do specific diagnoses determine higher levels of expenditure? Do people who make more money pay more for healthcare? If there is such a relationship, is it linear? Is the growth faster or slower than linear? Do certain diagnoses, while not individually affecting expenditure much, combine to have a large effect, perhaps larger than the sum of their individual pieces?

Looking more closely at the data will, I hope, help me validate or refute interesting hypotheses like these, and possibly generate many more. The first step is to try and visualize the possible relationships, which will be in the next post.

In the last post, I managed to construct the MEPS dataset as a csv file. Here's the result, for reference. (Warning: the file is large: about 250MB.)

The next thing to do is to take a look at the dataset, and see if we need to do any more processing. Since we have it written to disk as a csv file, I can import it using pandas' read_csv function. I may have to do a bit of plotting when I explore the dataset, so I'll go ahead and load all the packages I imagine I'll need for this.

Here's the boilerplate:

Let's load the dataset as a pandas dataframe and see what it contains:

First let's take a look at what features we have:

The first thing I get from this is that I need to learn what these features actually are, so I'll have to look up the MEPS survey methodology or documentation in order to figure out what these features actually are. There are 1792 features for each record, so this may be a bit involved.

The second thing I get out of looking at the column names is that they really do need a bit of cleanup. For example, the first feature is called "DWELLING UNIT ID", but the second is called " PERSON NUMBER" (note the extra space at the beginning). This seems to be a running problem in the column names. There's also a column right at the end with no name.

There are 1792 column names to fix, which is entirely too many to fix by hand.

On the other hand, it looks like the cleanup consists of a couple of simple rules:

Remove the initial space in the feature name
Pop the last entry from the column names

Let's write a function to do the first job:

Let's try this on our dataframe and see if it does the first task.

So far so good; now I just need to get rid of the last column.

Great, that worked. Alright, let's roll that into the clean_columns function:

Perfect.

One thing worth remarking on is that I've actually extracted the variable description rather than the variable name in this step. This was a deliberate decision, as I find it much more intuitive to work with variable names like "MOST RECENT COLONOSCOPY (>39) - R5/3" rather than the variable name "CLNTST53". Those of you familiar with the constraints of variable names in SQL databases will object that this is going to cause trouble in the future. I decided to go with the more descriptive variable names, and write a coding-decoding function to help move data in and out of SQL databases a bit later in the project. Overall, I found this to be a reasonable choice.

To this end, let me go ahead and get rid of the spaces in the variables names before that causes any bugs down the road.

Let's run the data columns through this function:

Great!

Next, I'd like to learn a bit more about what these features are. The main thing to keep in mind is that having a large number of features will make our model very unwieldy.

The crucial thing for this model is for it to be as useful as possible, and if I have to collect 1700 features or so about a person before I can give them any useful feedback, this model won't be used by anyone.

Secondly, for the purposes of my project, I'm most concerned with how the demographic, economic, or health status of a consumer drives cost, since these are the things that can inform how we choose to design the market for healthcare services.

To this end, let's take a look at what's available at the MEPS page that points to the dataset. This page gives us two helpful pieces of information: the first is the documentation for the dataset and the second is the codebook, which we've already used in constructing the csv file.

While the codebook contains the entire listing of variables in the dataset and the start-end points in the ASCII file, it does not tell us what values the variables actually take. To be fair, the variable names are fairly descriptive so we might be able to make educated guesses here, but a quick query of one of the columns in the dataset we've constructed:

raises the question: what do the values of this variable mean? It has a max of 2, and a min of -9. That doesn't correspond to any prior I currently have over the values of a diagnosis of high cholesterol. I'd expect to see things like positive, negative, not tested, and so on. Let's take a more detailed look:

The most frequently populated values are 2, 1 and -1, but there are also records with values -8, -7 and -9.

This warrants a closer look at the documentation. A quick search for "cholesterol" turns up:

High Cholesterol

Questions about high cholesterol were asked of persons aged 18 or older. Consequently, persons aged 17 or younger were coded as “Inapplicable” (-1) on these variables. CHOLDX ascertained whether the person had ever been diagnosed as having high cholesterol. Through 2007, a person-level variable (CHLAGE) indicated the age of diagnosis for high cholesterol on the Person-Level Use PUF. The age of diagnosis for high cholesterol (CHOLAGED) is included in this file. This variable is top-coded to 85 years of age.

which is somewhat helpful, but doesn't explain the use of the values -8 and -9.

Further digging into the documentation yields the following table:

2.2 Reserved Codes

The following reserved code values are used:

Value     Definition
-1 INAPPLICABLE     Question was not asked due to skip pattern
-2 DETERMINED IN PREVIOUS ROUND     Question was not asked in round because there was no change in current main job since previous round
-7 REFUSED     Question was asked and respondent refused to answer question
-8 DK     Question was asked and respondent did not know answer
-9 NOT ASCERTAINED     Interviewer did not record the data
-10 HOURLY WAGE >= $76.96     Hourly wage was top-coded for confidentiality
-13 INITIAL WAGE IMPUTED     Hourly wage was previously imputed so an updated wage is not included in this file

This is helpful; it indicates that the most relevant levels for our analysis ought to be the following three:

1. Positive diagnosis of high cholesterol
2. Negative diagnosis of high cholesterol
3. Inapplicable to respondent

and given the low population in levels below -1, we can either aggregate these data points with level==-1, or discard them entirely. In the following analysis, we choose to discard them.

A final question that needs to be asked at this step is which level corresponds to a positive diagnosis of high cholesterol? A close look at the documentation indicates that level==1 corresponds to having the condition, and level==2 corresponds to not having the condition. Knowing which level is which doesn't particularly affect the computations that generate our models going forward, but it is important for being able to interpret the models.

As a quick sanity check, let's get a visualization of the distribution of the data between these levels:

Of the respondents for whom the question was applicable, it looks like the ratio of level==1 to level==2 is 7/17 or about 2/5. This is close to the CDC estimates of the prevalence of high cholesterol amongst American adults: http://www.cdc.gov/dhdsp/data_statistics/fact_sheets/fs_cholesterol.htm .

Next, instead of using level==1 or level==2 to describe the values that a categorical feature takes, I'd like to translate the levels of these categorical variables into their descriptions so my plots are easier to read.

Here's a dictionary that implements this using the documentation quoted above:

Let's write a function that implements this as necessary for a feature.

The features whose levels I need to interpret in this way are all categorical features; the documentation indicates that the continuous variables aren't generally coded (apart from top-coding to preserve privacy). Accordingly, it would help to be able to automatically determine whether the feature I'm looking at is categorical or continuous. I'll try implementing the following rough rule of thumb: if a feature has less than 10-20 levels, I'll assume it's categorical, and otherwise I'll assume it's continuous.

Let's classify the features in the dataset:

Finally, I'll write a function that interprets the levels of the features in any categorical variable:

Here's an example of how it transforms the levels in a plot:

That visualization was a nice way of learning something about the data, and confirming a hunch we had about how to interpret it. In fact, simple visual plots of different features in the dataset are a crucial way of getting some initial insight into the data, and helping us prioritize the most fruitful approaches for modeling the data.

To fix ideas, let's scope down what we want to go into our model.

Explicitly, we'd like to be able to answer the following questions: which (if any) of the following features about a person determine what they will pay for healthcare in a year?

1. Demographic information:
   1. Age
   2. Sex
   3. Race/Ethnicity
   4. Marital Status
   5. Educational attainment
   6. Census region
2. Income information (these are all top-coded in the dataset for confidentiality reasons):
3. Idiosyncratic health status variables:
   1. High blood pressure, including multiple diagnoses
   2. Heart disease (including coronary heart disease, angina, myocardial infarction, and other unspecified heart disease)
   3. Stroke
   4. Emphysema
   5. Chronic Bronchitis
   6. High cholesterol, including the age of diagnosis
   7. Cancer
   8. Diabetes
   9. Joint pain
   10. Arthritis
   11. Asthma
   12. Attention Deficit Hyperactivity Disorder/Attention Deficit Disorder (ADHD/ADD)
4. Adult Body Mass Index
5. Health Insurance Variables
   1. Public Insurance Indicators
   2. Private Insurance Indicators
   3. Any Insurance Indicators
   4. HMO plan
6. Utilization, expenditure and source of payments information

It's fair to say that this is a lot of information to process. One thing that stands out here is that there are a great number of features already recorded for us in this dataset; this is in sharp contrast to some projects I've done using text and/or social media data, in which the hard part of the modeling process is feature construction in some appropriate way.

Since there is such a rich variety of features, I imagine most of the feature engineering in this project will amount to bucketing variables and levels in appropriate ways to optimize explained variation, rather than constructing new features from scratch.

To this end, let's make a list of the exact features we want to start working with.

1. AGE AS OF 12/31/13 (EDITED/IMPUTED)
2. SEX
3. RACE/ETHNICITY (EDITED/IMPUTED)
4. MARITAL STATUS-12/31/13 (EDITED/IMPUTED)
5. FAMILY'S TOTAL INCOME
6. ADULT BODY MASS INDEX (>17) - RD 5/3
7. HIGH BLOOD PRESSURE DIAG (>17)
8. CORONARY HRT DISEASE DIAG (>17)
9. AGE OF DIAGNOSIS-HIGH BLOOD PRESSURE
10. MULT DIAG HIGH BLOOD PRESS (>17)
11. AGE OF DIAGNOSIS-CORONARY HEART DISEASE
12. ANGINA DIAGNOSIS (>17)
13. AGE OF DIAGNOSIS-ANGINA
14. HEART ATTACK (MI) DIAG (>17)
15. AGE OF DIAGNOSIS-HEART ATTACK(MI)
16. OTHER HEART DISEASE DIAG (>17)
17. AGE OF DIAGNOSIS-OTHER HEART DISEASE
18. STROKE DIAGNOSIS (>17)
19. AGE OF DIAGNOSIS-STROKE
20. EMPHYSEMA DIAGNOSIS (>17)
21. AGE OF DIAGNOSIS-EMPHYSEMA
22. CHRONC BRONCHITS LAST 12 MTHS (>17)-R3/1
23. CHRONC BRONCHITS LAST 12 MTHS (>17)-R5/3
24. HIGH CHOLESTEROL DIAGNOSIS (>17)
25. AGE OF DIAGNOSIS-HIGH CHOLESTEROL
26. CANCER DIAGNOSIS (>17)
27. CANCER DIAGNOSED - BLADDER (>17)
28. CANCER DIAGNOSED - BREAST (>17)
29. CANCER DIAGNOSED - CERVIX (>17)
30. CANCER DIAGNOSED - COLON (>17)
31. CANCER DIAGNOSED - LUNG (>17)
32. CANCER DIAGNOSED - LYMPHOMA (>17)
33. CANCER DIAGNOSED - MELANOMA (>17)
34. CANCER DIAGNOSED - OTHER (>17)
35. CANCER DIAGNOSED - PROSTATE (>17)
36. CANCER DIAGNOSED - SKIN-NONMELANO (>17)
37. CANCER DIAGNOSED-SKIN-UNKNOWN TYPE (>17)
38. CANCER DIAGNOSED - UTERUS (>17)
39. DIABETES DIAGNOSIS (>17)
40. AGE OF DIAGNOSIS-DIABETES
41. JOINT PAIN LAST 12 MONTHS (>17) - RD 3/1
42. JOINT PAIN LAST 12 MONTHS (>17) - RD 5/3
43. ARTHRITIS DIAGNOSIS (>17)
44. AGE OF DIAGNOSIS-ARTHRITIS
45. ASTHMA DIAGNOSIS
46. AGE OF DIAGNOSIS-ASTHMA
47. DOES PERSON STILL HAVE ASTHMA-RD3/1
48. DOES PERSON STILL HAVE ASTHMA - RD 5/3
49. ASTHMA ATTACK LAST 12 MOS - RD3/1
50. USED ACUTE PRES INHALER LAST 3 MOS-RD5/3
51. USED>3ACUTE CN PRES INH LAST 3 MOS-RD5/3
52. EVER USED PREV DAILY ASTHMA MEDS -RD 5/3
53. NOW TAKE PREV DAILY ASTHMA MEDS - RD 5/3
54. HAVE PEAK FLOW METER AT HOME - RD 5/3
55. EVER USED PEAK FLOW METER - RD 5/3
56. ADHDADD DIAGNOSIS (5-17)
57. AGE OF DIAGNOSIS-ADHD/ADD
58. PREGNANT DURING REF PERIOD - RD 3/1
59. LIMITATION IN PHYSICAL FUNCTIONING-RD3/1
60. EMPLOYMENT STATUS RD 3/1
61. HAS MORE THAN ONE JOB RD 3/1 INT DATE
62. SELF-EMPLOYED AT RD 3/1 CMJ
63. CHOICE OF HEALTH PLANS AT RD 3/1 CMJ
64. INDUSTRY GROUP RD 3/1 CMJ
65. UNION STATUS AT RD 3/1 CMJ
66. OCCUPATION GROUP RD 3/1 CMJ
67. HEALTH INSUR HELD FROM RD 3/1 CMJ (ED)
68. HEALTH INSUR OFFERED BY RD 3/1 CMJ (ED)
69. EMPLOYER OFFERS HEALTH INS RD 3/1 CMJ
70. FULL YEAR INSURANCE COVERAGE STATUS 2013
71. ANYONE IN RU HAVE FSA - R3/1
72. # WKS/MON WOUT HLTH INS PRV YR-PN 18 ONL
73. PRIVATE INSURANCE ANY TIME IN R5/R3
74. PUBLIC INS ANY TIME IN R5/R3
75. INSURED ANY TIME IN R3/1
76. ANY TIME COVERAGE BY STATE INS - R3/1
77. # OFFICE-BASED PROVIDER VISITS 13
78. # OUTPATIENT DEPT PROVIDER VISITS 13
79. # EMERGENCY ROOM VISITS 13
80. # NIGHTS IN HOSP FOR DISCHARGES 2013
81. FINAL PERSON WEIGHT 2013
82. TOTAL OFFICE-BASED EXP 13
83. TOTAL OUTPATIENT PROVIDER EXP 13
84. TOT HOSP IP FACILITY + DR EXP 13
85. TOTAL ER FACILITY + DR EXP 13
86. CENSUS REGION AS OF 12/31/13
87. EDUCATION RECODE (EDITED)

It's going to be extremely interesting to see if these have explanatory power for the statistical variation in spending on healthcare!

Speaking of which, what are the variables we're trying to predict?

1. TOTAL OFFICE-BASED EXP 13
2. TOTAL OUTPATIENT PROVIDER EXP 13
3. TOT HOSP IP FACILITY + DR EXP 13
4. TOTAL ER FACILITY + DR EXP 13

The dataset breaks out healthcare expenses into these four location based categories: office, outpatient, inpatient and ER.

Another modeling choice that needs to be made is whether to try and predict total expenses, or predict each category of expense individually. In my Insight project, I opted to try and predict them individually, to see if I could shed some light on whether the service location affected prices charged. There's reason to think this may be the case: ER sticker shock is widely reported in the news, more so than office visits.

Finally, it's time to reflect a bit on what we've done here. We've isolated explanatory variables that we hope and expect should be predictive of how much a person is charged for healthcare services in four different categories. That's great, but there's a very large number of predictors: over 80. Since the ultimate goal is to build a usable interface to the model that people can use to get a sense of their expected outlays on healthcare, it makes sense to refine and reduce this list. But how?

One way to think about this is to consider what information people are likely to be comfortable sharing about themselves that could help to predict expenses. Another is to focus on what the most predictive variables are, and then ask users to populate these variables in order to get a prediction.

In my Insight project, I opted to balance both these concerns by ranking the explanatory variables in terms of statistical variation explained, cutting off at a reasonable threshold (somewhere in between 10 and 20 variables, corresponding to 10-20 questions asked of the user), and also culling variables that would be particularly burdensome to the user to have to report.

The next thing to think about is: now that we have explanatory variables and explained variables selected, how do we best go about building a model?

My next step is to examine the data feature by feature, to get a better sense of what it looks like and to pick up hunches for what may explain the variation. Along the way I'm sure I'll develop some more detailed hypotheses about healthcare spending, and discover that some of what I implicitly assume to be true about the healthcare market isn't actually so.

First, let's look at the explained variables:

This is interesting. It looks like upto half of respondents spend less than $200 for office-based services in a year. The maximum though, is concerning: over a hundred thousand dollars? Let's take a closer look at that.

There are 4 people with office based expenses of over 100k. I'd like to take a look to see if these records are bad data, or if they should be kept in the dataset. My prior is that they should be kept in the dataset (the MEPS survey has been running for over 15 years by people who I can only assume know what they're doing). But it doesn't hurt to take a closer look.

Now, I don't want to look at all the extra features in the dataset beyond those that I've decided I care about, so let's go about trimming the dataframe down to contain just these features.

As a matter of good practice, let's write the current dataset to a file, in case we break something and need to come back to this checkpoint.

I've manually written the explanatory variables to a text file located at "/home/soumya/research/insight/insight_project/modeling/code/dsm_exog_.txt". I'm going to read this into memory:

I'm going to make sure there are no duplicates in my list:

For future reference, let's write this dataset to disk so we don't have to do all this over again.

Let's check this worked:

Great! Now let's look at those four data points more closely. N.B: Note the change of assignment to the name "data".

Ok, if you look at the feature called "# OFFICE-BASED PROVIDER VISITS 13", it looks like these people had 168, 164, 51 and 128 visits each. While this certainly seems anomalous, it's consistent with the high values for amounts billed for office-based expenses.

At this point it's unclear whether these points should be left in the dataset, or thrown out. On the one hand, one would hope that most people don't need over 160 visits to an office-based provider in a year. On the other hand, it's a well-documented fact that a small handful of very ill people account for the vast majority of healthcare costs in the U.S. every year (this is well known to third-party payors). This may make it important to keep in these outliers in the final model.

I'd like to take a closer look at these points. What was the average charge per visit for these respondents?

Between $\$650$ and $\$2700$. Presumably these were specialists, and expensive.

If these amounts are not fraudulent, I can imagine that it would be important to leave them in while modeling, to help account for the risk of large healthcare bills. Since there are only four respondents with such high levels of expenditure, let's take a closer look to see if we can find out what drove their spending.

I'll take a look at these respondents one by one:

It looks like the outliers consist of respondents who are chronically unwell, and this may be driving heavy treatment regimens.

Given the evidence, I'm inclined to leave these respondents in the dataset for now. While they are a small proportion of the population, their cases may be indicative of the financial consequences of unmanaged illness. I'll have to make sure their inclusion is explicitly mentioned in my model writeup, since I am making a call on whether these points are corrupt data.

Next, let's look at the overall distribution of expenditure on office based services.

The graph above is a histogram of billed charges in the dataset. This graph is difficult to read because of the scale. I might try and break it in to a few pieces: one consisting of data points near zero (perhaps one with zero and one with strictly positive values), another consisting of data points that are positive but not too large (say less than 50k or so) and the final piece consisting of the extremely large positive values.

Let's try this:

I want to point out here that the dataset consists of weighted points; each point carries a weight (the column is "FINAL_PERSON_WEIGHT_2013") such that when the dataset is scaled out pointwise by these weights, it gives a representative sampling of the U.S. population. Since I do have these weights given to me, I've used these in my histograms to get a more accurate sense of the distribution of these features in the U.S. population.

This graph clearly shows the extent of the zero clustering. It looks like over 80m people are estimated to have expenses of less than $\$10$ - this is about one-fourth of the population.

Next, let's take a look at the same range, but without the mass at 0.

This shows that the modal nonzero expense for office services is about $\$100-\$200$ dollars, with about 32 million people falling into that bucket.

Let's look further out.

I've used a bin width of $\$100$ again, to keep things consistent. It looks like about 240m people spend less than $\$1000$ a year for office based services. Significant numbers of people (about 50m people) spend between $\$1000$ and $\$3000$ on office based services. After that, each $\$100$ bin contains relatively few people, but these are at fairly high levels of expenditure for an average person or family.

There is quite a bit of mass in the range $\$5000 - \$10,000$. These people, at 8.5 million people, or about $2.7\%$ of the population, are paying a lot for office based healthcare services. It would be interesting to know how much of expenditure at this level is discretionary or non-essential, and how much of it is for relatively routine healthcare - perhaps for badly managed chronic illnesses, etc. In terms of informing policy, there's nothing wrong with the former, but since chronic illness tends to affect a broad section of people - not all of whom can afford to pay upto $\$10,000$ a year for healthcare services, there may be something to worry about if a lot of this high expenditure is for non-discretionary services.

Let's look a bit further:

There is significant mass at these higher levels of expenditure - it's a small proportion of the population (4.5 million, or about $1.4\%$ of the population), but the severity of the costs is non-trivial. It will be very interesting to see if we can find out more about what characterizes people who tend to pay this much for healthcare in a year. Do they have terminal illnesses? Badly managed chronic illnesses? Severe accidents?

Let's look at the tail end of the distribution:

There are about 1.8 million people in this range, and the severity of the costs they face is astonishingly high. Again, it will be extremely interesting to know what about these people drives this level of expenditure.

Next, let's take a look at the other three predicted variables, in the following order:

1. Outpatient services
2. Inpatient services
3. Emergency Room services

Clearly, there's a lot of zero clustering going on here as well. Let's break it out into pieces like we did for office based services, and see what we find.

It looks like the vast majority of respondents pay less than $\$10$ for outpatient services.

Let's see what happens if we drop the zero payments:

Of the order of 300m people incur charges of a couple of hundred dollars for outpatient services in the year.

Let's zoom in a bit more to that area and see if we can figure out what the modal charge is.

The modal charge seems to be $\$30-\$40$ dollars for outpatient services.

Let's now look further to the right to see what the high expenditure cases look like:

There are about 5.6m people in this range, which is about $1.75\%$ of the population. There are about 1.1m people, or $0.34\%$ of respondents that are paying over $\$3,000$ in this range.

Let's take a look at the rest of the data.

Not many people here - about 77,000. The maximum amount paid by any respondent is $\$26,703$, so we won't find mass any further to the right.

It looks like outpatient expenditure follows a similar pattern to office based expenditure, although the severity is less extreme.

Let's turn to inpatient expenditure.

Alright, so we have a max of about $\$320,000$. That's extremely high!

Let's go through the same steps as before:

Again, most people pay less than $\$10$.

What if we drop the mass at 0?

Let's move out a bit further:

Ok, quite a lot of mass in this range. Charges in the thousands of dollars seem more likely in the context of inpatient services, since this usually includes procedures like surgeries.

Let's look out past $\$10,000$.

This is consistent with the kinds of amounts we're used to hearing about being charged for surgeries, etc. About 4m people seem to have inpatient services costing over $\$20,000$; these amounts are quite severely high. Let's see if it gets worse for anyone.

For the 1.4m people in this range, the amounts charged are more than the median household's annual income in the United States. Anyone that falls into this range - and it looks like many people do - will be at great financial risk.

I'd really like to know if the people who end up in this range have anything in common, or whether it's mostly random. Are these life-threatening accidents? Terminal illness? Badly managed chronic illness? Being able to tell the difference, if there is any, could be very informative for healthcare policy, and for advising healthcare consumers about their financial risk.

Let's move on to the final predicted variable, ER expenditure.

The maximum charge for ER services is $\$63,179$. This is a lot lower than the maximum for inpatient services, but still very high.

Let's look at this close up.

Fewer people at $\$0$ than for outpatient or inpatient expenditure, but more than for office based expenditure.

Let's drop the $\$0$ charges and see what it looks like.

There are more people in this range than there were for inpatient services, though less than for office based services and outpatient services. A side-by-side comparison:

While the overall utilization of ER services is lower than for office based services at each cost level, there seem to be many more people using ER services at the higher end of this range than either outpatient or inpatient services.

It's hard to say whether this means anything at this point, but it is interesting nonetheless.

Medical emergencies often come with significant sticker shock, so I'm curious to see what this looks like at much higher levels of expenditure.